

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <pthread.h>
#include <ctype.h>
#include "rand.h"
 
#define NUM 10
 
// Variaveis globais para os parametros de entrada do metódo, não são modificadas após lidas.
// TODO Trocar doubles para big numbers, não esqueca as alocações de memorias para os ponteiros  
double s,e,r,sigma,tempo;
// Variaveis globais para os resultados do método
double *trials, media = 0.0, stddev = 0.0;
// Numero de iterações (e threads)
int size,t_size;
 
//struct BoxMullerState state;
 
 
 
//Funcao que calcula media dos valores de um vetor
double calc_media() {
    int i;
    double result = 0.0;
    for (i=0; i < size; i++) {
        result = result + trials[i];
    }
   
    return result/size;
}
 
//Funcao que calcula o desvio padrao de um vetor
double calc_stddev() {
    int i;
    double result = 0.0, temp;
    for (i=0; i< size; i++) {
        temp = trials[i] - media;
        result += temp*temp;
    }  
    return sqrt(result/size);
}
 
//Algoritmo Black-Scholes com Monte Carlo
void *black_scholes() {
    double temp, **result;
    int i;
    struct BoxMullerState state;
    initBoxMullerState(&state);
    result = (double**) malloc(sizeof(double*)*t_size);
    for (i=0;i<t_size;i++) {
        result[i] = (double *) malloc(sizeof(double));
        temp = s*exp( (r-0.5*sigma*sigma)*tempo + sigma*sqrt(tempo)*boxMullerRandom(&state) );
        *result[i] = exp(-r*tempo)*fmax(temp-e,0);
    }
    return result;
}
 
int main(int argc, char *argv[]) {
 
    pthread_t *thread;
    pthread_attr_t attr;
    double confwidth, confmin, confmax, **thread_return;
    int i, j, x;
    double temp;
    int num = NUM;


    if(argc > 1)
        num = atoi(argv[1]);


    //initBoxMullerState(&state);  
    //srand(0); // Inicializa RNG, talvez trocar 0 por funçao time()
   
    //Leitura da entrada
    scanf("%lf",&s);
    scanf("%lf",&e);
    scanf("%lf",&r);
    scanf("%lf",&sigma);
    scanf("%lf",&tempo);
    scanf("%d",&size);
 
    //nthread = size;
    //nthread = ceil(nthread/NUM);
    temp = size;
    t_size = ceil(temp/num);
 
    //Aloca vetor de trials e thread IDs
    trials = (double*) malloc(sizeof(double)*(int) size);
    thread = (pthread_t*) malloc(sizeof(pthread_t)*num);
 
    pthread_attr_init(&attr);
    pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
     
    //Cria as threas que executarao as iteracoes do black scholes
    for (i=0; i < num; i++) {    
       pthread_create(&thread[i], &attr, &black_scholes, NULL);
    }
 
    pthread_attr_destroy(&attr);
 
    for (i=0; i< num; i++) {
        x = pthread_join(thread[i], (void**) &thread_return);
        if (x != 0) {
            printf("Erro no join %d %d\n", x, i);
            exit(0);
        }
 
        for(j = 0; j < t_size; j++) {
            if (i*t_size+j < size) {
                trials[i*t_size+j] = *thread_return[j];
                free(thread_return[j]);
            }
        }
        free(thread_return);
    //printf("Thread %d completo\n", i);
    }
 
    //Calcula media dos trials
    media = calc_media();
        //printf("Media: %lf\n",media);  
     
    //Desvio padrão dos trials
    stddev = calc_stddev();
    //printf("stddev: %lf\n",stddev);
   
    //Calcula intervalo de confianca  
    confwidth = 1.96*stddev/sqrt(size);
    confmin = media - confwidth;
    confmax = media + confwidth;
 
    // Saida formatada do programa
    printf("S        %lf\nE        %lf\nr        %lf\nsigma    %lf\nT        %lf\nM        %d\n",s,e,r,sigma,tempo,size);
    printf("Confidence interval: (%lf,%lf)\n",confmin,confmax);
 
    free(trials);
    free(thread);
    return 0;            
}

